SGE Behavior Project

# Load behavioral data
load('/Users/noah/Desktop/SGE_data_objects/sge_phase12.Rdata')
data$Fq9_SGE2 <- data$Fq9_SGE2[c(1:5012,5144:5791),]
data$Pp10_SGE2 <- data$Pp10_SGE2[c(1:662,797:7334),]

## load functions & intros object 
load('/Users/noah/Desktop/SGE_data_objects/behav_functions.RData')
load('/Users/noah/Desktop/SGE_data_objects/datafornoah.RData')

## number of focals per group
for (x in names(data)) intros[which(intros[,"group"]==x),"num.obs"]=length(levels(as.factor((data[[x]][,"obs.no"]))))

# pull out social behaviors and dominance interactions

social<-lapply(data,convert.soc)
elo<-lapply(data,convert.elo)
groom<-lapply(social,function(x){x[which(x[,5]%in%c("GM","G-","ZZ")),]})
dominance<-lapply(data,convert.dom)
prox<-lapply(social,function(x){x[which(x[,5]%in%c("PX","LB","ZZ")),]})

intros$finalElo <- NA
for (i in 1:90){
  id=tolower(intros[i,"ID"])
  group=intros[i,"group"]
  intros[i,"finalElo"]=extract_elo(elo[[group]],ID=id)
}

intros<-intros[order(intros$phase,intros$group),]

## Calculate SGE I hierarchy metrics and generate a dataframe containing all metrics
sge.I.hierarchy.metrics <- matrix(nrow = 18, ncol = 6)
count <- 1
for (i in unique(names(dominance))){
  tmp <- as.data.frame(dominance[i])
  tmp <- tmp[order(tmp[,1]),]
  tmp <- subset(tmp, tmp[,3] != tmp[,4])
  seqcheck(tmp[,3],tmp[,4],tmp[,1])
  elo<-elo.seq(tmp[,3],tmp[,4],tmp[,1])
  elo.mat <- creatematrix(elo)
  a <- extract_elo(elo)
  b <- steeptest(elo.mat, rep = 1000)$Stp
  c <- devries(elo.mat)$`h-modified`
  d <- dci(elo.mat)
  e <- phi(elo.mat)
  f <- ttri(elo.mat)$ttri
  g <- stab_elo(elo)
  sge.I.hierarchy.metrics[count,] <- c(b,c,d,e,f,g)
  count <- count + 1
}
sge.I.hierarchy.metrics <- as.data.frame(sge.I.hierarchy.metrics)
sge.I.hierarchy.metrics$phase <- rep(c(1,2), times = c(9,9))
rownames(sge.I.hierarchy.metrics) <- names(dominance)
colnames(sge.I.hierarchy.metrics) <- c('steepness','linearity','DC','phi','tt','stability','phase')
sge.I.hierarchy.metrics$group <- rownames(sge.I.hierarchy.metrics)

Supplementary information about hierarchies and behavioral variables

Descriptive statistics of hierarchies

Here I’m calculating a variety of metrics that describe the nature of the dominance hierarchies in SGE I. These include steepness, linearity, direction concordance, phi, triangle transitivity, and stability (descriptions of these hierarchy metrics can be found here).

Correlations between hierarchy metrics.
Panel 1 (both phase); Panel 2 (phase 1); panel 3 (phase 2). Triangle transitivity is not included in the plots as all groups had a value of 1. Similarly, all groups but one had a value of 1 for linearity.

# plot correlations between elos, DS, and ordinal ranks for SGE I
ggpairs(as.data.frame(sge.I.hierarchy.metrics[,c(1,3:4,6)]), title = 'Phases 1 + 2')

ggpairs(as.data.frame(sge.I.hierarchy.metrics[1:9,c(1,3:4,6)]), title = 'Phase 1')

ggpairs(as.data.frame(sge.I.hierarchy.metrics[10:18,c(1,3:4,6)]), title = 'Phase 2')

Distribution of hierarchy metrics across groups. Phi is plotted separately as it’s max value is 0.5 rather than 1 as is the case with the other metrics.

metric_melt <- sge.I.hierarchy.metrics[,c(1:3,5:6)]
metric_melt <- reshape2::melt(metric_melt)
phi_melt <- sge.I.hierarchy.metrics[,4]
phi_melt <- reshape2::melt(phi_melt)
p1 <- ggplot(metric_melt, aes(x=as.factor(variable), y = value)) + 
  geom_boxplot() +
  xlab('hierarchy metric')
p2 <- ggplot(phi_melt, aes(y = value)) + 
  geom_boxplot() +
  xlab('phi')
plot_grid(p1,p2)

Behavioral Variables

Here I’m going to show the correlations between all the behavioral variables being used in the study.

# read in behavioral data
sge.I.behavioral.v2 <- read.delim("~/Desktop/SGE_data_objects/sge.I.behavioral.v2.txt")
sge.I.behavioral.v2$rowOrder <- 1:nrow(sge.I.behavioral.v2)
sge.I.behavioral.v2$aggDiff <- sge.I.behavioral.v2$raw.aggGiven.rate - sge.I.behavioral.v2$raw.aggRec.rate

# PC1.agg and PC1.groom
agg_vector <- c('raw.aggRec.contact.rate', 'raw.aggRec.nocontact.rate', 'raw.aggGiven.contact.rate', 'raw.aggGiven.nocontact.rate', 'aggDiff')
tmp1 <- subset(sge.I.behavioral.v2, phase == '1') 
tmp2 <- subset(sge.I.behavioral.v2, phase == '2')
pca1 <- prcomp(as.data.frame(scale(tmp1[,agg_vector])), scale = T)
pca2 <- prcomp(as.data.frame(scale(tmp2[,agg_vector])), scale = T)
sge.I.behavioral.v2$PC1.agg <- NA
sge.I.behavioral.v2[which(sge.I.behavioral.v2$phase=='1'),'PC1.agg'] <- pca1$x[,1]
sge.I.behavioral.v2[which(sge.I.behavioral.v2$phase=='2'),'PC1.agg'] <- pca2$x[,1]

groom_vector <- c('raw.groomRec.rate', 'raw.groomGiven.rate')
tmp1 <- subset(sge.I.behavioral.v2, phase == '1') 
tmp2 <- subset(sge.I.behavioral.v2, phase == '2')
pca1 <- prcomp(as.data.frame(scale(tmp1[,groom_vector])), scale = T)
pca2 <- prcomp(as.data.frame(scale(tmp2[,groom_vector])), scale = T)
sge.I.behavioral.v2$PC1.groom <- NA
sge.I.behavioral.v2[which(sge.I.behavioral.v2$phase=='1'),'PC1.groom'] <- pca1$x[,1]
sge.I.behavioral.v2[which(sge.I.behavioral.v2$phase=='2'),'PC1.groom'] <- pca2$x[,1]

# group-center and scale behavioral variables
sge.I.behavioral.v2$overall.agg.received <- sge.I.behavioral.v2$raw.aggRec.contact.rate + sge.I.behavioral.v2$raw.aggRec.nocontact.rate
sge.I.behavioral.v2$overall.agg.given <- sge.I.behavioral.v2$raw.aggGiven.contact.rate + sge.I.behavioral.v2$raw.aggGiven.nocontact.rate
sge.I.behavioral.v2$overall.groom <- sge.I.behavioral.v2$raw.groomRec.rate + sge.I.behavioral.v2$raw.groomGiven.rate

sge.I.behavioral.v2 <- sge.I.behavioral.v2[order(sge.I.behavioral.v2$group),]

group.center <- function(behaveVar){
  test <- vector()
  for (i in unique(sge.I.behavioral.v2$group)){
    tmp <- subset(sge.I.behavioral.v2, group == i)
    tmp2 <- scale(tmp[,behaveVar], scale=F, center=T)
    test <- rbind(test, tmp2)
  }
  sge.I.behavioral.v2[,paste('gc.',behaveVar,sep='')] <- as.numeric(scale(test))
  return(sge.I.behavioral.v2)
}

behaviorList <- c('overall.groom',
                  'overall.agg.given',
                  'overall.agg.received',
                  'raw.aggRec.contact.rate',
                  'raw.aggRec.nocontact.rate',
                  'raw.aggGiven.contact.rate',
                  'raw.aggGiven.nocontact.rate',
                  'raw.groomRec.rate',
                  'raw.groomGiven.rate',
                  'aggDiff')

for (i in behaviorList){
  sge.I.behavioral.v2 <- group.center(i)
}

sge.I.behavioral.v2 <- sge.I.behavioral.v2[order(sge.I.behavioral.v2$rowOrder),]

Elo vs. agg variables

p1 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.aggDiff, group = phase)) +
  geom_point(aes(color=as.factor(phase))) +
  xlab('rank (elo)') +
  ylab('aggDiff') +
  stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
  ggtitle('') +
  scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
  theme_minimal() 
p2 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=PC1.agg, group = phase)) +
  geom_point(aes(color=as.factor(phase))) +
  xlab('rank (elo)') +
  ylab('PC1.agg') +
  stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
  ggtitle('') +
  scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
  theme_minimal()
p3 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.overall.agg.given, group = phase)) +
  geom_point(aes(color=as.factor(phase))) +
  xlab('rank (elo)') +
  ylab('Overall agg given') +
  stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
  ggtitle('') +
  scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
  theme_minimal()
p4 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.overall.agg.received, group = phase)) +
  geom_point(aes(color=as.factor(phase))) +
  xlab('rank (elo)') +
  ylab('Overall agg received') +
  stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
  ggtitle('') +
  scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
  theme_minimal()
plot_grid(p1,p2,p3,p4, ncol = 2)

Elo vs. groom variables

p1 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.raw.groomGiven.rate, group = phase)) +
  geom_point(aes(color=as.factor(phase))) +
  xlab('rank (elo)') +
  ylab('Overall groom given') +
  stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
  ggtitle('') +
  scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
  theme_minimal()
p2 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.raw.groomRec.rate, group = phase)) +
  geom_point(aes(color=as.factor(phase))) +
  xlab('rank (elo)') +
  ylab('Overall groom received') +
  stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
  ggtitle('') +
  scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
  theme_minimal()
plot_grid(p1,p2, ncol = 2)

All pairwise comparisons

ggpairs(sge.I.behavioral.v2[,c(5,26:27,31:40)],
        upper = list(continuous = wrap("cor", size = 2)),
        lower = list(continuous = wrap("points", size = 0.5))) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Now we can look at the PVE and loadings for PC1.agg and PC1.groom.
First we’ll look at PC1.agg

## Generate PCA info for PC1.agg variables for SGE I and II both phases

sgeIbehavior <- read.table('/Users/noah/Desktop/SGE_data_objects/sge.I.behavioral.v2.txt', header = T)
sgeIbehavior$aggDiff <- sgeIbehavior$raw.aggGiven.rate - sgeIbehavior$raw.aggRec.rate
sgeIIbehavior <- read.table('/Users/noah/Desktop/SGE_data_objects/sge.II.behavioral.v2.txt', header = T)
sgeIIbehavior$aggDiff <- sgeIIbehavior$raw.aggGiven.rate - sgeIIbehavior$raw.aggRec.rate

agg_vector <- c('raw.aggRec.contact.rate', 'raw.aggRec.nocontact.rate', 'raw.aggGiven.contact.rate', 'raw.aggGiven.nocontact.rate', 'aggDiff')

tmp1 <- subset(sgeIbehavior, phase == '1') 
tmp2 <- subset(sgeIbehavior, phase == '2') 
tmp3 <- subset(sgeIIbehavior, phase == '1')
tmp4 <- subset(sgeIIbehavior, phase == '2')

pca1 <- prcomp(as.data.frame(scale(tmp1[,agg_vector])), scale = T)
pca2 <- prcomp(as.data.frame(scale(tmp2[,agg_vector])), scale = T)
pca3 <- prcomp(as.data.frame(scale(tmp3[,agg_vector])), scale = T)
pca4 <- prcomp(as.data.frame(scale(tmp4[,agg_vector])), scale = T)

pov_table <- matrix(nrow = 20, ncol = 2)
pov_table[,1] <- rep(c('I_I','I_II','II_I','II_II'), times=c(5,5,5,5))
pov_table[1:5,2] <- summary(pca1)$importance[2,]
pov_table[6:10,2] <- summary(pca2)$importance[2,]
pov_table[11:15,2] <- summary(pca3)$importance[2,]
pov_table[16:20,2] <- summary(pca4)$importance[2,]
pov_table <- as.data.frame(pov_table)
pov_table$PC <- rep(1:5, times=4)

pov_table <- pov_table[1:10,]
pov_table$V1 <- as.factor(rep(c(1,2), times=c(5,5)))
ggplot(pov_table, aes(x=PC, y = as.numeric(V2), color = V1)) +
  geom_line(size=2) +
  geom_point(size=2) +
  xlab('PC') +
  ylab('Proportion of variance explained by PC') +
  labs(color = "Phase") +
  theme_minimal() +
  scale_color_viridis(discrete = T, option = 'E')

I_I <- cbind(reshape2::melt(pca1$rotation), 'I_I')
I_II <- cbind(reshape2::melt(pca2$rotation), 'I_II')
II_I <- cbind(reshape2::melt(pca3$rotation), 'II_I')
II_II <- cbind(reshape2::melt(pca4$rotation), 'II_II')
loading_list <- list(I_I,I_II,II_I,II_II)
loading_table <- data.table::rbindlist(loading_list, use.names = F)
colnames(loading_table) <- c('behavior', 'PC', 'loading', 'exp_phase')

loading_table <- loading_table[1:50,]
loading_table$exp_phase <- as.factor(rep(c(1,2), times = c(25,25)))
ggplot(loading_table, aes(PC, abs(loading))) +
  geom_linerange(aes(x = PC, ymin = 0, ymax = abs(loading), group = behavior), color = "lightgray", size = 1.5, position = position_dodge(0.3)) +
  geom_point(aes(color = behavior),position = position_dodge(0.3), size = 3) +
  facet_grid(rows = vars(exp_phase)) +
  theme_minimal() +
  scale_color_viridis(discrete = T, option = 'E')

Now we can look at PC1.groom

## Generate PCA info for PC1.agg variables for SGE I and II both phases

sgeIbehavior <- read.table('/Users/noah/Desktop/SGE_data_objects/sge.I.behavioral.v2.txt', header = T)
sgeIbehavior$aggDiff <- sgeIbehavior$raw.aggGiven.rate - sgeIbehavior$raw.aggRec.rate
sgeIIbehavior <- read.table('/Users/noah/Desktop/SGE_data_objects/sge.II.behavioral.v2.txt', header = T)
sgeIIbehavior$aggDiff <- sgeIIbehavior$raw.aggGiven.rate - sgeIIbehavior$raw.aggRec.rate

groom_vector <- c('raw.groomRec.rate', 'raw.groomGiven.rate')

tmp1 <- subset(sgeIbehavior, phase == '1') 
tmp2 <- subset(sgeIbehavior, phase == '2') 
tmp3 <- subset(sgeIIbehavior, phase == '1')
tmp4 <- subset(sgeIIbehavior, phase == '2')

pca1 <- prcomp(as.data.frame(scale(tmp1[,groom_vector])), scale = T)
pca2 <- prcomp(as.data.frame(scale(tmp2[,groom_vector])), scale = T)
pca3 <- prcomp(as.data.frame(scale(tmp3[,groom_vector])), scale = T)
pca4 <- prcomp(as.data.frame(scale(tmp4[,groom_vector])), scale = T)

pov_table <- matrix(nrow = 8, ncol = 2)
pov_table[,1] <- rep(c('I_I','I_II','II_I','II_II'), times=c(2,2,2,2))
pov_table[1:2,2] <- summary(pca1)$importance[2,]
pov_table[3:4,2] <- summary(pca2)$importance[2,]
pov_table[5:6,2] <- summary(pca3)$importance[2,]
pov_table[7:8,2] <- summary(pca4)$importance[2,]
pov_table <- as.data.frame(pov_table)
pov_table$PC <- rep(1:2, times=4)

pov_table <- pov_table[1:4,]
pov_table$V1 <- as.factor(rep(c(1,2), times = c(2,2)))
ggplot(pov_table, aes(x=PC, y = as.numeric(V2), color = V1)) +
  geom_line(size=2) +
  geom_point(size=2) +
  xlab('PC') +
  ylab('Proportion of variance explained by PC') +
  labs(color = "Phase") +
  theme_minimal() +
  scale_color_viridis(discrete = T, option = 'E')

I_I <- cbind(reshape2::melt(pca1$rotation), 'I_I')
I_II <- cbind(reshape2::melt(pca2$rotation), 'I_II')
II_I <- cbind(reshape2::melt(pca3$rotation), 'II_I')
II_II <- cbind(reshape2::melt(pca4$rotation), 'II_II')
loading_list <- list(I_I,I_II,II_I,II_II)
loading_table <- data.table::rbindlist(loading_list, use.names = F)
colnames(loading_table) <- c('behavior', 'PC', 'loading', 'exp_phase')

loading_table <- loading_table[1:8,]
loading_table$exp_phase <- as.factor(rep(c(1,2), times = c(4,4)))
ggplot(loading_table, aes(PC, abs(loading))) +
  geom_linerange(aes(x = PC, ymin = 0, ymax = abs(loading), group = behavior), color = "lightgray", size = 1.5, position = position_dodge(0.3)) +
  geom_point(aes(color = behavior),position = position_dodge(0.3), size = 3) +
  facet_grid(rows = vars(exp_phase)) +
  theme_minimal() +
  scale_color_viridis(discrete = T, option = 'E')

As you can see the loadings are not meaningful with only two variables in the PCA, but the first PC explains ~60% of the variance. Not sure is PC1.groom is a valid composite variable.

Assessment of data quality

One metric we will use to assess the data quality is to split the data and see how correlated the behavioral variables are between the two subsets, where we’d expect the two subsets to be well correlated if we’re capturing ‘true’ behavioral rates.

Because the subsets need to maintain the row structure in order to calculate grooming times, the data were split by obs, alternating across the study duration so both subsets contain data from across the entire study.

##caclulate grooming bouts
## Here I'm going to 1) split the data in half, 2) regenerate groom objects, 3) calculate grooming rates for each half, and 4) see how correlated they are
## Can't split by taking alternating rows becasue the row stricture needs to be maintained in order to calculate mins of grooming, so need to split by something that's unrelated to grooming, I'll sort alphabetically, split, and the resort to original order in each subset (by row order)

data_row <- lapply(names(data),function(x){data[[x]]$rowOrder <- 1:nrow(data[[x]]); return(data[[x]])})
names(data_row) <- names(data)

data_half1 <- lapply(names(data_row),function(x){evens <- seq(2,max(unique(data[[x]]$obs.no)),2);return(data_row[[x]][which(data_row[[x]][,3] %in% evens),])})
names(data_half1) <- names(data_row)

data_half2 <- lapply(names(data_row),function(x){odds <- seq(1,max(unique(data[[x]]$obs.no)),2);return(data_row[[x]][which(data_row[[x]][,3] %in% odds),])})
names(data_half2) <- names(data_row)

social_1 <- lapply(data_half1,convert.soc)
groom_1 <- lapply(social_1,function(x){x[which(x[,5]%in%c("GM","G-","ZZ")),]})
names(groom_1) <- names(data)
social_2 <- lapply(data_half2,convert.soc)
groom_2 <- lapply(social_2,function(x){x[which(x[,5]%in%c("GM","G-","ZZ")),]})
names(groom_2) <- names(data)

groom.function <- function(groomDF){
groom.tot1<-data.frame(matrix(ncol=19,nrow=360))
colnames(groom.tot1)<-c("ID1","ID2","group","no.obs","bouts","ID1.groom","ID2.groom","min","ID1.groom.min","ID2.groom.min","ID1.rank","ID2.rank","proximity","contact.min","CSI","elo.ID1","elo.ID2","ID1.agg","ID2.agg")
groom.tot1[,"ID1"]<-rep(intros[,"ID"],4)
for (i in 1:180) {groom.tot1[i,"group"]<-intros[which(intros[1:45,"ID"]==groom.tot1[i,1]),"group"]}; rm(i)
for (i in 181:360) {groom.tot1[i,"group"]<-intros[46:90,][which(intros[46:90,"ID"]==groom.tot1[i,1]),"group"]}; rm(i)
groom.tot1<-groom.tot1[order(groom.tot1[,3],groom.tot1[,1]),];rownames(groom.tot1)<-NULL

for (i in unique(intros[,4])) {intros[1:45,][which(intros[1:45,"ID"]==i),"group"]->grp; 
  groom.tot1[grep("SGE1",groom.tot1$group),][which(groom.tot1[grep("SGE1",groom.tot1$group),"ID1"]==i),2]<-setdiff(intros[1:45,][which(intros[1:45,"group"]==grp),"ID"],i)}

for (i in unique(intros[,4])) {intros[46:90,][which(intros[46:90,"ID"]==i),"group"]->grp; 
  groom.tot1[grep("SGE2",groom.tot1$group),][which(groom.tot1[grep("SGE2",groom.tot1$group),"ID1"]==i),2]<-setdiff(intros[46:90,][which(intros[46:90,"group"]==grp),"ID"],i)}


groom.tot<-groom.tot1; rm(groom.tot1)
groom.tot[,1]<-as.character(groom.tot[,1]); groom.tot[,2]<-as.character(groom.tot[,2]); groom.tot[,3]<-as.character(groom.tot[,3]); rm(grp); rm(i)

##caclulate grooming bouts
groom.tot$ID1.intro.no<-NA; groom.tot$ID2.intro.no<-NA
for (id in 1:nrow(groom.tot)){
  id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0;grp<-groom.tot[id,"group"]
  b1=sum(g[,5]=="GM"&g[,3]==id1 & g[,4]==id2)
  b2=sum(g[,5]=="GM"&g[,4]==id1 & g[,3]==id2)
  groom.tot[id,5]=b1+b2;groom.tot[id,6]=b1;groom.tot[id,7]=b2
  groom.tot[id,4]=intros[which(intros[,"ID"]==id2&intros$group==grp),"num.obs"]
  groom.tot[id,"ID1.rank"]=intros[which(intros[,"ID"]==id1&intros$group==grp),"finalElo"]
  groom.tot[id,"ID2.rank"]=intros[which(intros[,"ID"]==id2&intros$group==grp),"finalElo"]
  groom.tot[id,"ID1.intro.no"]=intros[which(intros[,"ID"]==id1&intros$group==grp),"intro.no"]
  groom.tot[id,"ID2.intro.no"]=intros[which(intros[,"ID"]==id2&intros$group==grp),"intro.no"]  
}; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)

groom.tot$ID1.intro.date<-NA; groom.tot$ID2.intro.date<-NA
for (id in 1:nrow(groom.tot)){
  id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0;grp<-groom.tot[id,"group"]
  groom.tot[id,"ID1.intro.date"]=intros[which(intros[,"ID"]==id1&intros$group==grp),"intro.date"]
  groom.tot[id,"ID2.intro.date"]=intros[which(intros[,"ID"]==id2&intros$group==grp),"intro.date"]  
}; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)

## calculate who ended the grooming bout
groom.tot$ID1.end.bout<-NA; groom.tot$ID2.end.bout<-NA
for (id in 1:nrow(groom.tot)){
  id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0
  b1=sum(g[,5]=="G-"&g[,3]==id1 & g[,4]==id2)
  b2=sum(g[,5]=="G-"&g[,4]==id1 & g[,3]==id2)
  groom.tot[id,"ID1.end.bout"]=b1
  groom.tot[id,"ID2.end.bout"]=b2
}; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)

groom.tot$ID1.approach<-NA
groom.tot$ID2.approach<-NA
##caclulate approaches
for (id in 1:nrow(groom.tot)){
  id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-prox[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0
  b1=sum(g[,5]=="PX"&g[,3]==id1 & g[,4]==id2)
  b2=sum(g[,5]=="PX"&g[,4]==id1 & g[,3]==id2)
  groom.tot[id,"ID1.approach"]=b1;groom.tot[id,"ID2.approach"]=b2
}; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)

## calculate grooming times
groom.tot$ID1.groom.min<-NA
groom.tot$ID2.groom.min<-NA
groom.tot$min<-NA

for (i in 1:nrow(groom.tot)){
  id1<-groom.tot[i,"ID1"]; id2<-groom.tot[i,"ID2"]; t1=0;t2=0
  g<-groomDF[[as.character(groom.tot[i,"group"])]]
  if (id1%in%c(g[,3],g[,4])&id2%in%c(g[,3],g[,4])){
    g<-g[which((g[,4]==id1 & g[,3]==id2 )|(g[,3]==id1 & g[,4]==id2)|g[,3]=="ZZ"),]
    for (r in 1:nrow(g)){
      if (g[r,5]=="GM" & g[r,3]==id1 & g[r,4]==id2) t1=t1+(g[r+1,2]-g[r,2])
      if (g[r,5]=="GM" & g[r,4]==id1 & g[r,3]==id2) t2=t2+(g[r+1,2]-g[r,2])
    }
    groom.tot[i,"ID1.groom.min"]=t1
    groom.tot[i,"ID2.groom.min"]=t2
    groom.tot[i,"min"]=t1+t2}
}
rm(t1);rm(t2);rm(i);rm(id1);rm(id2);rm(r);rm(g)

###########
phase <- groom.tot$group
phase <- substr(phase,(nchar(phase)+1)-4,nchar(phase))
groom.tot$phase <- phase

## generate a data frame with minutes of grooming (total), elo, # of obs, and group, per individual
## phase I
tmp.minutes <- vector()
tmp.minutes.given <- vector()
tmp.minutes.received <- vector()
tmp.names <- vector()
tmp.ranks <- vector()
tmp.obs <- vector()
tmp.group <- vector()
for (i in unique(groom.tot$ID1)){
  tmp <- subset(groom.tot, groom.tot$phase == 'SGE1')
  tmp2 <- sum(subset(tmp, tmp$ID1 == i)$min)
  tmp3 <- sum(subset(tmp, tmp$ID1 == i)$ID1.groom.min)
  tmp4 <- sum(subset(tmp, tmp$ID1 == i)$ID2.groom.min)
  tmp.minutes <- cbind(tmp.minutes,tmp2)
  tmp.minutes.given <- cbind(tmp.minutes.given,tmp3)
  tmp.minutes.received <- cbind(tmp.minutes.received,tmp4)
  tmp.names <- cbind(tmp.names,i)
  tmp.ranks <- cbind(tmp.ranks,subset(tmp, tmp$ID1 == i)$ID1.rank[1])
  tmp.obs <- cbind(tmp.obs,subset(tmp, tmp$ID1 == i)$no.obs[1])
  tmp.group <- cbind(tmp.group,subset(tmp, tmp$ID1 == i)$group[1])
}  
sge.I.groom.rate.df <- as.data.frame(t(rbind(tmp.names,tmp.ranks,tmp.minutes,tmp.minutes.given,tmp.minutes.received,tmp.obs,tmp.group)))
names(sge.I.groom.rate.df) <- c('ID','elo','min','min.given','min.received','no.obs','group')
sge.I.groom.rate.df$elo <- as.numeric(as.character(sge.I.groom.rate.df$elo))
sge.I.groom.rate.df$min <- as.numeric(as.character(sge.I.groom.rate.df$min))
sge.I.groom.rate.df$min.given <- as.numeric(as.character(sge.I.groom.rate.df$min.given))
sge.I.groom.rate.df$min.received <- as.numeric(as.character(sge.I.groom.rate.df$min.received))
sge.I.groom.rate.df$no.obs <- as.numeric(as.character(sge.I.groom.rate.df$no.obs))
sge.I.groom.rate.df$phase <- 1

## now for phase II
tmp.minutes <- vector()
tmp.minutes.given <- vector()
tmp.minutes.received <- vector()
tmp.names <- vector()
tmp.ranks <- vector()
tmp.obs <- vector()
tmp.group <- vector()
for (i in unique(groom.tot$ID1)){
  tmp <- subset(groom.tot, groom.tot$phase == 'SGE2')
  tmp2 <- sum(subset(tmp, tmp$ID1 == i)$min)
  tmp3 <- sum(subset(tmp, tmp$ID1 == i)$ID1.groom.min)
  tmp4 <- sum(subset(tmp, tmp$ID1 == i)$ID2.groom.min)
  tmp.minutes <- cbind(tmp.minutes,tmp2)
  tmp.minutes.given <- cbind(tmp.minutes.given,tmp3)
  tmp.minutes.received <- cbind(tmp.minutes.received,tmp4)
  tmp.names <- cbind(tmp.names,i)
  tmp.ranks <- cbind(tmp.ranks,subset(tmp, tmp$ID1 == i)$ID1.rank[1])
  tmp.obs <- cbind(tmp.obs,subset(tmp, tmp$ID1 == i)$no.obs[1])
  tmp.group <- cbind(tmp.group,subset(tmp, tmp$ID1 == i)$group[1])
}  
sge.II.groom.rate.df <- as.data.frame(t(rbind(tmp.names,tmp.ranks,tmp.minutes,tmp.minutes.given,tmp.minutes.received,tmp.obs,tmp.group)))
names(sge.II.groom.rate.df) <- c('ID','elo','min','min.given','min.received','no.obs','group')
sge.II.groom.rate.df$elo <- as.numeric(as.character(sge.II.groom.rate.df$elo))
sge.II.groom.rate.df$min <- as.numeric(as.character(sge.II.groom.rate.df$min))
sge.II.groom.rate.df$min.given <- as.numeric(as.character(sge.II.groom.rate.df$min.given))
sge.II.groom.rate.df$min.received <- as.numeric(as.character(sge.II.groom.rate.df$min.received))
sge.II.groom.rate.df$no.obs <- as.numeric(as.character(sge.II.groom.rate.df$no.obs))
sge.II.groom.rate.df$phase <- 2

sge.groom.rate.df <- rbind(sge.I.groom.rate.df,sge.II.groom.rate.df)
return(sge.groom.rate.df)
}

groom.df.1 <- groom.function(groom_1)
groom.df.2 <- groom.function(groom_2)

# turn minutes into rates
groom.df.1$min.rate <- (groom.df.1$min / groom.df.1$no.obs) * 2
groom.df.1$min.given.rate <- (groom.df.1$min.given / groom.df.1$no.obs) * 2
groom.df.1$min.received.rate <- (groom.df.1$min.received / groom.df.1$no.obs) * 2
groom.df.2$min.rate <- (groom.df.2$min / groom.df.2$no.obs) * 2
groom.df.2$min.given.rate <- (groom.df.2$min.given / groom.df.2$no.obs) * 2
groom.df.2$min.received.rate <- (groom.df.2$min.received / groom.df.2$no.obs) * 2

tmp <- cor.test(groom.df.1$min.rate,groom.df.2$min.rate)
p1 <- ggplot(groom.df.1, aes(x=min.rate,y=groom.df.2$min.rate)) +
  geom_point(alpha = 0.6, color=brocolors("crayons")["Bittersweet"]) +
  #geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
  stat_smooth(method='lm', color=brocolors("crayons")["Bittersweet"]) +
  xlab('grooming rate subset 1') +
  ylab('grooming rate subset 2') +
  ggtitle('total grooming rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) #+
  #stat_poly_eq(formula = y~x, 
  #             aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
  #             parse = TRUE)

tmp <- cor.test(groom.df.1$min.given.rate,groom.df.2$min.given.rate)
p2 <- ggplot(groom.df.1, aes(x=min.given.rate,y=groom.df.2$min.given.rate)) +
  geom_point(alpha = 0.6, color=brocolors("crayons")["Bittersweet"]) +
  #geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
  stat_smooth(method='lm', color=brocolors("crayons")["Bittersweet"]) +
  xlab('grooming rate subset 1') +
  ylab('grooming rate subset 2') +
  ggtitle('grooming given rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) #+
  #stat_poly_eq(formula = y~x, 
  #             aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
  #             parse = TRUE)

tmp <- cor.test(groom.df.1$min.received.rate,groom.df.2$min.received.rate)
p3 <- ggplot(groom.df.1, aes(x=min.received.rate,y=groom.df.2$min.received.rate)) +
  geom_point(alpha = 0.6, color=brocolors("crayons")["Bittersweet"]) +
  #geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
  stat_smooth(method='lm', color=brocolors("crayons")["Bittersweet"]) +
  xlab('grooming rate subset 1') +
  ylab('grooming rate subset 2') +
  ggtitle('grooming received rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) #+
  #stat_poly_eq(formula = y~x, 
  #             aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
  #             parse = TRUE)
p1

p2

p3

Ok, looks pretty strongly correlated, now we’ll look at the agg behaviors

SGEI.intros <- intros

# aggression received
agg.function <- function(socialDF){
  sge.groom.rate.df <- groom.df.1
  aggs <- prox<-lapply(socialDF,function(x){x[which(x[,5]%in%c("TN","TC","AT","CH")),]})
  aggs.contact <- lapply(socialDF,function(x){x[which(x[,5]%in%c("TC","AT")),]})
  aggs.nocontact <- lapply(socialDF,function(x){x[which(x[,5]%in%c("TN","CH")),]})
  agg.received<-lapply(aggs,function(x){aggregate(x[,4],by=list(x[,4]),length)})
  agg.given<-lapply(aggs,function(x){aggregate(x[,3],by=list(x[,3]),length)})
  agg.contact.received<-lapply(aggs.contact,function(x){aggregate(x[,4],by=list(x[,4]),length)})
  agg.contact.given<-lapply(aggs.contact,function(x){aggregate(x[,3],by=list(x[,3]),length)})
  agg.nocontact.received<-lapply(aggs.nocontact,function(x){aggregate(x[,4],by=list(x[,4]),length)})
  agg.nocontact.given<-lapply(aggs.nocontact,function(x){aggregate(x[,3],by=list(x[,3]),length)})
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
  tmp3 <- agg.received[[i]]$Group.1
  tmp4 <- agg.received[[i]]$x
  tmp5 <- rep(i,length(tmp3))
  agg <- c(agg,tmp4)
  name <- c(name,tmp3)
  group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')

intro_sub <- SGEI.intros
intro_sub$agg.received <- 0
intro_sub$id.group <- paste(intro_sub$ID,intro_sub$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub$id.group)
intro_sub <- subset(intro_sub, intro_sub$id.group %in% tmp$id.group)
intro_sub <- intro_sub[match(tmp$id.group, intro_sub$id.group),]
intro_sub <- na.omit(intro_sub)
intro_sub$agg.received <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub$id.group))
tmp$agg.received <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub <- rbind.data.frame(intro_sub,tmp)
####

agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
  tmp3 <- agg.contact.received[[i]]$Group.1
  tmp4 <- agg.contact.received[[i]]$x
  tmp5 <- rep(i,length(tmp3))
  agg <- c(agg,tmp4)
  name <- c(name,tmp3)
  group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')

intro_sub2 <- SGEI.intros
intro_sub2$agg.contact.received <- 0
intro_sub2$id.group <- paste(intro_sub2$ID,intro_sub2$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub2$id.group)
intro_sub2 <- subset(intro_sub2, intro_sub2$id.group %in% tmp$id.group)
intro_sub2 <- intro_sub2[match(tmp$id.group, intro_sub2$id.group),]
intro_sub2$agg.contact.received <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub2$id.group))
tmp$agg.contact.received <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub2 <- rbind.data.frame(intro_sub2,tmp)

####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
  tmp3 <- agg.nocontact.received[[i]]$Group.1
  tmp4 <- agg.nocontact.received[[i]]$x
  tmp5 <- rep(i,length(tmp3))
  agg <- c(agg,tmp4)
  name <- c(name,tmp3)
  group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
####
intro_sub3 <- SGEI.intros
intro_sub3$agg.nocontact.received <- 0
intro_sub3$id.group <- paste(intro_sub3$ID,intro_sub3$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub3$id.group)
intro_sub3 <- subset(intro_sub3, intro_sub3$id.group %in% tmp$id.group)
intro_sub3 <- intro_sub3[match(tmp$id.group, intro_sub3$id.group),]
intro_sub3$agg.nocontact.received <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub3$id.group))
tmp$agg.nocontact.received <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub3 <- rbind.data.frame(intro_sub3,tmp)

intro_sub <- intro_sub[order(intro_sub$id.group),]
intro_sub2 <- intro_sub2[order(intro_sub2$id.group),]
intro_sub3 <- intro_sub3[order(intro_sub3$id.group),]

sge.I.behavioral.df2 <- cbind(intro_sub[,c('ID','group','phase','agg.received','num.obs')],intro_sub2$agg.contact.received,intro_sub3$agg.nocontact.received)

sge.I.behavioral.df2 <- sge.I.behavioral.df2[order(sge.I.behavioral.df2$phase,sge.I.behavioral.df2$group,sge.I.behavioral.df2$ID),]
sge.groom.rate.df <- sge.groom.rate.df[order(sge.groom.rate.df$phase,sge.groom.rate.df$group,sge.groom.rate.df$ID),]

tmp.final.behav.df <- sge.I.behavioral.df2[,c(4,6:7)]


# aggression given
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
  tmp3 <- agg.given[[i]]$Group.1
  tmp4 <- agg.given[[i]]$x
  tmp5 <- rep(i,length(tmp3))
  agg <- c(agg,tmp4)
  name <- c(name,tmp3)
  group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')

intro_sub <- SGEI.intros
intro_sub$agg.given <- 0
intro_sub$id.group <- paste(intro_sub$ID,intro_sub$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub$id.group)
intro_sub <- subset(intro_sub, intro_sub$id.group %in% tmp$id.group)
intro_sub <- intro_sub[match(tmp$id.group, intro_sub$id.group),]
intro_sub$agg.given <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub$id.group))
tmp$agg.given <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub <- rbind.data.frame(intro_sub,tmp)
####

agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
  tmp3 <- agg.contact.given[[i]]$Group.1
  tmp4 <- agg.contact.given[[i]]$x
  tmp5 <- rep(i,length(tmp3))
  agg <- c(agg,tmp4)
  name <- c(name,tmp3)
  group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')

intro_sub2 <- SGEI.intros
intro_sub2$agg.contact.given <- 0
intro_sub2$id.group <- paste(intro_sub2$ID,intro_sub2$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub2$id.group)
intro_sub2 <- subset(intro_sub2, intro_sub2$id.group %in% tmp$id.group)
intro_sub2 <- intro_sub2[match(tmp$id.group, intro_sub2$id.group),]
intro_sub2$agg.contact.given <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub2$id.group))
tmp$agg.contact.given <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub2 <- rbind.data.frame(intro_sub2,tmp)

####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
  tmp3 <- agg.nocontact.given[[i]]$Group.1
  tmp4 <- agg.nocontact.given[[i]]$x
  tmp5 <- rep(i,length(tmp3))
  agg <- c(agg,tmp4)
  name <- c(name,tmp3)
  group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
####
intro_sub3 <- SGEI.intros
intro_sub3$agg.nocontact.given <- 0
intro_sub3$id.group <- paste(intro_sub3$ID,intro_sub3$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub3$id.group)
intro_sub3 <- subset(intro_sub3, intro_sub3$id.group %in% tmp$id.group)
intro_sub3 <- intro_sub3[match(tmp$id.group, intro_sub3$id.group),]
intro_sub3$agg.nocontact.given <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub3$id.group))
tmp$agg.nocontact.given <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub3 <- rbind.data.frame(intro_sub3,tmp)

intro_sub <- intro_sub[order(intro_sub$id.group),]
intro_sub2 <- intro_sub2[order(intro_sub2$id.group),]
intro_sub3 <- intro_sub3[order(intro_sub3$id.group),]

sge.I.behavioral.df2 <- cbind(intro_sub[,c('ID','group','phase','agg.given','num.obs')],intro_sub2$agg.contact.given,intro_sub3$agg.nocontact.given)


sge.I.behavioral.df2 <- sge.I.behavioral.df2[order(sge.I.behavioral.df2$phase,sge.I.behavioral.df2$group,sge.I.behavioral.df2$ID),]
sge.groom.rate.df <- sge.groom.rate.df[order(sge.groom.rate.df$phase,sge.groom.rate.df$group,sge.groom.rate.df$ID),]

final.behav.df <- cbind(sge.I.behavioral.df2,tmp.final.behav.df)
colnames(final.behav.df) <- c('ID',
                              'group',
                              'phase',
                              'agg.given',
                              'num.obs',
                              'agg.contact.given',
                              'agg.nocontact.given',
                              'agg.received',
                              'agg.contact.received',
                              'agg.nocontact.received')

final.behav.df$agg.given.rate <- (final.behav.df$agg.given / final.behav.df$num.obs) * 2
final.behav.df$agg.contact.given.rate <- (final.behav.df$agg.contact.given / final.behav.df$num.obs) * 2
final.behav.df$agg.nocontact.given.rate <- (final.behav.df$agg.nocontact.given / final.behav.df$num.obs) * 2
final.behav.df$agg.received.rate <- (final.behav.df$agg.received / final.behav.df$num.obs) * 2
final.behav.df$agg.contact.received.rate <- (final.behav.df$agg.contact.received / final.behav.df$num.obs) * 2
final.behav.df$agg.nocontact.received.rate <- (final.behav.df$agg.nocontact.received / final.behav.df$num.obs) * 2
return(final.behav.df)
}

agg.df.1 <- agg.function(social_1)
agg.df.2 <- agg.function(social_2)

tmp <- cor.test(agg.df.1$agg.given.rate,agg.df.2$agg.given.rate)
p1 <- ggplot(agg.df.1, aes(x=agg.given.rate,y=agg.df.2$agg.given.rate)) +
  geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
  #geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
  stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
  xlab('agg given rate subset 1') +
  ylab('agg given rate subset 2') +
  ggtitle('agg given rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) 

tmp <- cor.test(agg.df.1$agg.contact.given.rate,agg.df.2$agg.contact.given.rate)
p2 <- ggplot(agg.df.1, aes(x=agg.contact.given.rate,y=agg.df.2$agg.contact.given.rate)) +
  geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
  #geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
  stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
  xlab('agg contact given rate subset 1') +
  ylab('agg contact given rate subset 2') +
  ggtitle('agg contact given rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) 

tmp <- cor.test(agg.df.1$agg.nocontact.given.rate,agg.df.2$agg.nocontact.given.rate)
p3 <- ggplot(agg.df.1, aes(x=agg.nocontact.given.rate,y=agg.df.2$agg.nocontact.given.rate)) +
  geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
  #geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
  stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
  xlab('agg noncontact given rate subset 1') +
  ylab('agg noncontact given rate subset 2') +
  ggtitle('agg noncontact given rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) 

tmp <- cor.test(agg.df.1$agg.received.rate,agg.df.2$agg.received.rate)
p4 <- ggplot(agg.df.1, aes(x=agg.received.rate,y=agg.df.2$agg.received.rate)) +
  geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
  #geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
  stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
  xlab('agg received rate subset 1') +
  ylab('agg received rate subset 2') +
  ggtitle('agg received rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) 

tmp <- cor.test(agg.df.1$agg.contact.received.rate,agg.df.2$agg.contact.received.rate)
p5 <- ggplot(agg.df.1, aes(x=agg.contact.received.rate,y=agg.df.2$agg.contact.received.rate)) +
  geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
  #geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
  stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
  xlab('agg contact received rate subset 1') +
  ylab('agg contact received rate subset 2') +
  ggtitle('agg contact received rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) 

tmp <- cor.test(agg.df.1$agg.nocontact.received.rate,agg.df.2$agg.nocontact.received.rate)
p6 <- ggplot(agg.df.1, aes(x=agg.nocontact.received.rate,y=agg.df.2$agg.nocontact.received.rate)) +
  geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
  #geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
  stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
  xlab('agg noncontact received rate subset 1') +
  ylab('agg noncontact received rate subset 2') +
  ggtitle('agg noncontact received rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) 

p1

p2

p3

p4

p5

p6

Everything looks pretty well correlated. Perhaps as expected the noncontact agg measures are better correlated than contact agg measures.

Now let’s look at how the grooming and agg behaviors are distributed across a sliding window of obs. Here I’ve calculated grooming and agg data for 29 windows representing ~6 months of data collection and ~ 60 hours of data per group. The window size is 3 obs.

agg.function <- function(socialDF){
  sge.groom.rate.df <- groom.df.1
  SGEI.intros_tmp <- intros_tmp
  
  aggs <- prox<-lapply(socialDF,function(x){x[which(x[,5]%in%c("TN","TC","AT","CH")),]})
  aggs.contact <- lapply(socialDF,function(x){x[which(x[,5]%in%c("TC","AT")),]})
  
  #aggs.contact[[2]] <- aggs.contact[[1]]
  #aggs.contact[[5]] <- aggs.contact[[1]]
  #aggs.contact[[16]] <- aggs.contact[[1]]
  
  for (i in 1:18){
  if (nrow(aggs.contact[[i]]) == 0){
    aggs.contact[[i]][1,1:6] <- c('2013-10-10',1,1,1,1,1)
  }
  }
  
  
  for (i in 1:18){
    if (nrow(aggs[[i]]) == 0){
      aggs[[i]][1,1:6] <- c('2013-10-10',1,1,1,1,1)
    }
  }
  aggs.nocontact <- lapply(socialDF,function(x){x[which(x[,5]%in%c("TN","CH")),]})
  
  for (i in 1:18){
    if (nrow(aggs.nocontact[[i]]) == 0){
      aggs.nocontact[[i]][1,1:6] <- c('2013-10-10',1,1,1,1,1)
    }
  }
  
  agg.received<-lapply(aggs,function(x){aggregate(x[,4],by=list(x[,4]),length)})
  agg.given<-lapply(aggs,function(x){aggregate(x[,3],by=list(x[,3]),length)})
  agg.contact.received<-lapply(aggs.contact,function(x){aggregate(x[,4],by=list(x[,4]),length)})
  agg.contact.given<-lapply(aggs.contact,function(x){aggregate(x[,3],by=list(x[,3]),length)})
  agg.nocontact.received<-lapply(aggs.nocontact,function(x){aggregate(x[,4],by=list(x[,4]),length)})
  agg.nocontact.given<-lapply(aggs.nocontact,function(x){aggregate(x[,3],by=list(x[,3]),length)})
  agg <- vector()
  name <- vector()
  group <- vector()
  for (i in unique(SGEI.intros_tmp$group)){
    tmp3 <- agg.received[[i]]$Group.1
    tmp4 <- agg.received[[i]]$x
    tmp5 <- rep(i,length(tmp3))
    agg <- c(agg,tmp4)
    name <- c(name,tmp3)
    group <- c(group,tmp5)
  }
  tmp <- as.data.frame(cbind(agg,name,group))
  tmp <- tmp[!(tmp$name=="ZZ"),]
  tmp$agg <- as.numeric(as.character(tmp$agg))
  tmp$name <- toupper(tmp$name)
  tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
  
  intro_sub <- SGEI.intros_tmp
  intro_sub$agg.received <- 0
  intro_sub$id.group <- paste(intro_sub$ID,intro_sub$group, sep='_')
  tmp <- subset(tmp, tmp$id.group %in% intro_sub$id.group)
  intro_sub <- subset(intro_sub, intro_sub$id.group %in% tmp$id.group)
  intro_sub <- intro_sub[match(tmp$id.group, intro_sub$id.group),]
  intro_sub <- na.omit(intro_sub)
  intro_sub$agg.received <- tmp$agg
  # add zero rows back in
  SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
  tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub$id.group))
  tmp$agg.received <- 0
  #tmp <- tmp[,c(1:8,10,9)]
  intro_sub <- rbind.data.frame(intro_sub,tmp)
  ####
  
  agg <- vector()
  name <- vector()
  group <- vector()
  for (i in unique(SGEI.intros_tmp$group)){
    tmp3 <- agg.contact.received[[i]]$Group.1
    tmp4 <- agg.contact.received[[i]]$x
    tmp5 <- rep(i,length(tmp3))
    agg <- c(agg,tmp4)
    name <- c(name,tmp3)
    group <- c(group,tmp5)
  }
  tmp <- as.data.frame(cbind(agg,name,group))
  tmp <- tmp[!(tmp$name=="ZZ"),]
  tmp$agg <- as.numeric(as.character(tmp$agg))
  tmp$name <- toupper(tmp$name)
  tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
  
  intro_sub2 <- SGEI.intros_tmp
  intro_sub2$agg.contact.received <- 0
  intro_sub2$id.group <- paste(intro_sub2$ID,intro_sub2$group, sep='_')
  tmp <- subset(tmp, tmp$id.group %in% intro_sub2$id.group)
  intro_sub2 <- subset(intro_sub2, intro_sub2$id.group %in% tmp$id.group)
  intro_sub2 <- intro_sub2[match(tmp$id.group, intro_sub2$id.group),]
  intro_sub2$agg.contact.received <- tmp$agg
  # add zero rows back in
  SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
  tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub2$id.group))
  tmp$agg.contact.received <- 0
  #tmp <- tmp[,c(1:8,10,9)]
  intro_sub2 <- rbind.data.frame(intro_sub2,tmp)
  
  ####
  agg <- vector()
  name <- vector()
  group <- vector()
  for (i in unique(SGEI.intros_tmp$group)){
    tmp3 <- agg.nocontact.received[[i]]$Group.1
    tmp4 <- agg.nocontact.received[[i]]$x
    tmp5 <- rep(i,length(tmp3))
    agg <- c(agg,tmp4)
    name <- c(name,tmp3)
    group <- c(group,tmp5)
  }
  tmp <- as.data.frame(cbind(agg,name,group))
  tmp <- tmp[!(tmp$name=="ZZ"),]
  tmp$agg <- as.numeric(as.character(tmp$agg))
  tmp$name <- toupper(tmp$name)
  tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
  ####
  intro_sub3 <- SGEI.intros_tmp
  intro_sub3$agg.nocontact.received <- 0
  intro_sub3$id.group <- paste(intro_sub3$ID,intro_sub3$group, sep='_')
  tmp <- subset(tmp, tmp$id.group %in% intro_sub3$id.group)
  intro_sub3 <- subset(intro_sub3, intro_sub3$id.group %in% tmp$id.group)
  intro_sub3 <- intro_sub3[match(tmp$id.group, intro_sub3$id.group),]
  intro_sub3$agg.nocontact.received <- tmp$agg
  # add zero rows back in
  SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
  tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub3$id.group))
  tmp$agg.nocontact.received <- 0
  #tmp <- tmp[,c(1:8,10,9)]
  intro_sub3 <- rbind.data.frame(intro_sub3,tmp)
  
  intro_sub <- intro_sub[order(intro_sub$id.group),]
  intro_sub2 <- intro_sub2[order(intro_sub2$id.group),]
  intro_sub3 <- intro_sub3[order(intro_sub3$id.group),]
  
  sge.I.behavioral.df2 <- cbind(intro_sub[,c('ID','group','phase','agg.received','num.obs')],intro_sub2$agg.contact.received,intro_sub3$agg.nocontact.received)
  
  sge.I.behavioral.df2 <- sge.I.behavioral.df2[order(sge.I.behavioral.df2$phase,sge.I.behavioral.df2$group,sge.I.behavioral.df2$ID),]
  sge.groom.rate.df <- sge.groom.rate.df[order(sge.groom.rate.df$phase,sge.groom.rate.df$group,sge.groom.rate.df$ID),]
  
  tmp.final.behav.df <- sge.I.behavioral.df2[,c(4,6:7)]
  
  
  # aggression given
  agg <- vector()
  name <- vector()
  group <- vector()
  for (i in unique(SGEI.intros_tmp$group)){
    tmp3 <- agg.given[[i]]$Group.1
    tmp4 <- agg.given[[i]]$x
    tmp5 <- rep(i,length(tmp3))
    agg <- c(agg,tmp4)
    name <- c(name,tmp3)
    group <- c(group,tmp5)
  }
  tmp <- as.data.frame(cbind(agg,name,group))
  tmp <- tmp[!(tmp$name=="ZZ"),]
  tmp$agg <- as.numeric(as.character(tmp$agg))
  tmp$name <- toupper(tmp$name)
  tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
  
  intro_sub <- SGEI.intros_tmp
  intro_sub$agg.given <- 0
  intro_sub$id.group <- paste(intro_sub$ID,intro_sub$group, sep='_')
  tmp <- subset(tmp, tmp$id.group %in% intro_sub$id.group)
  intro_sub <- subset(intro_sub, intro_sub$id.group %in% tmp$id.group)
  intro_sub <- intro_sub[match(tmp$id.group, intro_sub$id.group),]
  intro_sub$agg.given <- tmp$agg
  # add zero rows back in
  SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
  tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub$id.group))
  tmp$agg.given <- 0
  #tmp <- tmp[,c(1:8,10,9)]
  intro_sub <- rbind.data.frame(intro_sub,tmp)
  ####
  
  agg <- vector()
  name <- vector()
  group <- vector()
  for (i in unique(SGEI.intros_tmp$group)){
    tmp3 <- agg.contact.given[[i]]$Group.1
    tmp4 <- agg.contact.given[[i]]$x
    tmp5 <- rep(i,length(tmp3))
    agg <- c(agg,tmp4)
    name <- c(name,tmp3)
    group <- c(group,tmp5)
  }
  tmp <- as.data.frame(cbind(agg,name,group))
  tmp <- tmp[!(tmp$name=="ZZ"),]
  tmp$agg <- as.numeric(as.character(tmp$agg))
  tmp$name <- toupper(tmp$name)
  tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
  
  intro_sub2 <- SGEI.intros_tmp
  intro_sub2$agg.contact.given <- 0
  intro_sub2$id.group <- paste(intro_sub2$ID,intro_sub2$group, sep='_')
  tmp <- subset(tmp, tmp$id.group %in% intro_sub2$id.group)
  intro_sub2 <- subset(intro_sub2, intro_sub2$id.group %in% tmp$id.group)
  intro_sub2 <- intro_sub2[match(tmp$id.group, intro_sub2$id.group),]
  intro_sub2$agg.contact.given <- tmp$agg
  # add zero rows back in
  SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
  tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub2$id.group))
  tmp$agg.contact.given <- 0
  #tmp <- tmp[,c(1:8,10,9)]
  intro_sub2 <- rbind.data.frame(intro_sub2,tmp)
  
  ####
  agg <- vector()
  name <- vector()
  group <- vector()
  for (i in unique(SGEI.intros_tmp$group)){
    tmp3 <- agg.nocontact.given[[i]]$Group.1
    tmp4 <- agg.nocontact.given[[i]]$x
    tmp5 <- rep(i,length(tmp3))
    agg <- c(agg,tmp4)
    name <- c(name,tmp3)
    group <- c(group,tmp5)
  }
  tmp <- as.data.frame(cbind(agg,name,group))
  tmp <- tmp[!(tmp$name=="ZZ"),]
  tmp$agg <- as.numeric(as.character(tmp$agg))
  tmp$name <- toupper(tmp$name)
  tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
  ####
  intro_sub3 <- SGEI.intros_tmp
  intro_sub3$agg.nocontact.given <- 0
  intro_sub3$id.group <- paste(intro_sub3$ID,intro_sub3$group, sep='_')
  tmp <- subset(tmp, tmp$id.group %in% intro_sub3$id.group)
  intro_sub3 <- subset(intro_sub3, intro_sub3$id.group %in% tmp$id.group)
  intro_sub3 <- intro_sub3[match(tmp$id.group, intro_sub3$id.group),]
  intro_sub3$agg.nocontact.given <- tmp$agg
  # add zero rows back in
  SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
  tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub3$id.group))
  tmp$agg.nocontact.given <- 0
  #tmp <- tmp[,c(1:8,10,9)]
  intro_sub3 <- rbind.data.frame(intro_sub3,tmp)
  
  intro_sub <- intro_sub[order(intro_sub$id.group),]
  intro_sub2 <- intro_sub2[order(intro_sub2$id.group),]
  intro_sub3 <- intro_sub3[order(intro_sub3$id.group),]
  
  sge.I.behavioral.df2 <- cbind(intro_sub[,c('ID','group','phase','agg.given','num.obs')],intro_sub2$agg.contact.given,intro_sub3$agg.nocontact.given)
  
  
  sge.I.behavioral.df2 <- sge.I.behavioral.df2[order(sge.I.behavioral.df2$phase,sge.I.behavioral.df2$group,sge.I.behavioral.df2$ID),]
  sge.groom.rate.df <- sge.groom.rate.df[order(sge.groom.rate.df$phase,sge.groom.rate.df$group,sge.groom.rate.df$ID),]
  
  final.behav.df <- cbind(sge.I.behavioral.df2,tmp.final.behav.df)
  colnames(final.behav.df) <- c('ID',
                                'group',
                                'phase',
                                'agg.given',
                                'num.obs',
                                'agg.contact.given',
                                'agg.nocontact.given',
                                'agg.received',
                                'agg.contact.received',
                                'agg.nocontact.received')
  
  final.behav.df$agg.given.rate <- (final.behav.df$agg.given / final.behav.df$num.obs) * 2
  final.behav.df$agg.contact.given.rate <- (final.behav.df$agg.contact.given / final.behav.df$num.obs) * 2
  final.behav.df$agg.nocontact.given.rate <- (final.behav.df$agg.nocontact.given / final.behav.df$num.obs) * 2
  final.behav.df$agg.received.rate <- (final.behav.df$agg.received / final.behav.df$num.obs) * 2
  final.behav.df$agg.contact.received.rate <- (final.behav.df$agg.contact.received / final.behav.df$num.obs) * 2
  final.behav.df$agg.nocontact.received.rate <- (final.behav.df$agg.nocontact.received / final.behav.df$num.obs) * 2
  return(final.behav.df)
}

#sliding window
windows.matrix <- matrix(nrow = 27, ncol = 3)
for (i in 1:27){
  windows.matrix[i,] <- seq(i,i+2,1)
}
windows.matrix <- cbind(windows.matrix, paste(windows.matrix[,1],':',windows.matrix[,3],sep=''))

for (i in 1:27){
  tmp <- lapply(names(data_row),function(x){
    tmpRange <- windows.matrix[i,4]
    tmp2 <- subset(data_row[[x]], obs.no %in% seq(i,i+2,1))
    return(tmp2)
  })
  names(tmp) <- names(data_row)
  assign(value = tmp, x = paste('data_window_',i,sep=''))
}  

groom.function2 <- function(groomDF){
  groom.tot1<-data.frame(matrix(ncol=19,nrow=360))
  colnames(groom.tot1)<-c("ID1","ID2","group","no.obs","bouts","ID1.groom","ID2.groom","min","ID1.groom.min","ID2.groom.min","ID1.rank","ID2.rank","proximity","contact.min","CSI","elo.ID1","elo.ID2","ID1.agg","ID2.agg")
  groom.tot1[,"ID1"]<-rep(intros_tmp[,"ID"],4)
  for (i in 1:180) {groom.tot1[i,"group"]<-intros_tmp[which(intros_tmp[1:45,"ID"]==groom.tot1[i,1]),"group"]}; rm(i)
  for (i in 181:360) {groom.tot1[i,"group"]<-intros_tmp[46:90,][which(intros_tmp[46:90,"ID"]==groom.tot1[i,1]),"group"]}; rm(i)
  groom.tot1<-groom.tot1[order(groom.tot1[,3],groom.tot1[,1]),];rownames(groom.tot1)<-NULL
  
  for (i in unique(intros_tmp[,4])) {intros_tmp[1:45,][which(intros_tmp[1:45,"ID"]==i),"group"]->grp; 
    groom.tot1[grep("SGE1",groom.tot1$group),][which(groom.tot1[grep("SGE1",groom.tot1$group),"ID1"]==i),2]<-setdiff(intros_tmp[1:45,][which(intros_tmp[1:45,"group"]==grp),"ID"],i)}
  
  for (i in unique(intros_tmp[,4])) {intros_tmp[46:90,][which(intros_tmp[46:90,"ID"]==i),"group"]->grp; 
    groom.tot1[grep("SGE2",groom.tot1$group),][which(groom.tot1[grep("SGE2",groom.tot1$group),"ID1"]==i),2]<-setdiff(intros_tmp[46:90,][which(intros_tmp[46:90,"group"]==grp),"ID"],i)}
  
  
  groom.tot<-groom.tot1; rm(groom.tot1)
  groom.tot[,1]<-as.character(groom.tot[,1]); groom.tot[,2]<-as.character(groom.tot[,2]); groom.tot[,3]<-as.character(groom.tot[,3]); rm(grp); rm(i)
  
  ##caclulate grooming bouts
  groom.tot$ID1.intro.no<-NA; groom.tot$ID2.intro.no<-NA
  for (id in 1:nrow(groom.tot)){
    id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0;grp<-groom.tot[id,"group"]
    b1=sum(g[,5]=="GM"&g[,3]==id1 & g[,4]==id2)
    b2=sum(g[,5]=="GM"&g[,4]==id1 & g[,3]==id2)
    groom.tot[id,5]=b1+b2;groom.tot[id,6]=b1;groom.tot[id,7]=b2
    groom.tot[id,4]=intros_tmp[which(intros_tmp[,"ID"]==id2&intros_tmp$group==grp),"num.obs"]
    groom.tot[id,"ID1.rank"]=intros_tmp[which(intros_tmp[,"ID"]==id1&intros_tmp$group==grp),"finalElo"]
    groom.tot[id,"ID2.rank"]=intros_tmp[which(intros_tmp[,"ID"]==id2&intros_tmp$group==grp),"finalElo"]
    groom.tot[id,"ID1.intro.no"]=intros_tmp[which(intros_tmp[,"ID"]==id1&intros_tmp$group==grp),"intro.no"]
    groom.tot[id,"ID2.intro.no"]=intros_tmp[which(intros_tmp[,"ID"]==id2&intros_tmp$group==grp),"intro.no"]  
  }#; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
  
  groom.tot$ID1.intro.date<-NA; groom.tot$ID2.intro.date<-NA
  for (id in 1:nrow(groom.tot)){
    id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0;grp<-groom.tot[id,"group"]
    groom.tot[id,"ID1.intro.date"]=intros_tmp[which(intros_tmp[,"ID"]==id1&intros_tmp$group==grp),"intro.date"]
    groom.tot[id,"ID2.intro.date"]=intros_tmp[which(intros_tmp[,"ID"]==id2&intros_tmp$group==grp),"intro.date"]  
  }#; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
  
  ## calculate who ended the grooming bout
  groom.tot$ID1.end.bout<-NA; groom.tot$ID2.end.bout<-NA
  for (id in 1:nrow(groom.tot)){
    id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0
    b1=sum(g[,5]=="G-"&g[,3]==id1 & g[,4]==id2)
    b2=sum(g[,5]=="G-"&g[,4]==id1 & g[,3]==id2)
    groom.tot[id,"ID1.end.bout"]=b1
    groom.tot[id,"ID2.end.bout"]=b2
  }#; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
  
  groom.tot$ID1.approach<-NA
  groom.tot$ID2.approach<-NA
  ##caclulate approaches
  for (id in 1:nrow(groom.tot)){
    id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-prox[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0
    b1=sum(g[,5]=="PX"&g[,3]==id1 & g[,4]==id2)
    b2=sum(g[,5]=="PX"&g[,4]==id1 & g[,3]==id2)
    groom.tot[id,"ID1.approach"]=b1;groom.tot[id,"ID2.approach"]=b2
  }#; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
  
  ## calculate grooming times
  groom.tot$ID1.groom.min<-NA
  groom.tot$ID2.groom.min<-NA
  groom.tot$min<-NA
  
  for (i in 1:nrow(groom.tot)){
    id1<-groom.tot[i,"ID1"]; id2<-groom.tot[i,"ID2"]; t1=0;t2=0
    g<-groomDF[[as.character(groom.tot[i,"group"])]]
    if (id1%in%c(g[,3],g[,4])&id2%in%c(g[,3],g[,4])){
      g<-g[which((g[,4]==id1 & g[,3]==id2 )|(g[,3]==id1 & g[,4]==id2)|g[,3]=="ZZ"),]
      for (r in 1:nrow(g)){
        if (g[r,5]=="GM" & g[r,3]==id1 & g[r,4]==id2) t1=t1+(g[r+1,2]-g[r,2])
        if (g[r,5]=="GM" & g[r,4]==id1 & g[r,3]==id2) t2=t2+(g[r+1,2]-g[r,2])
      }
      groom.tot[i,"ID1.groom.min"]=t1
      groom.tot[i,"ID2.groom.min"]=t2
      groom.tot[i,"min"]=t1+t2}
  }
  #rm(t1);rm(t2);rm(i);rm(id1);rm(id2);rm(r);rm(g)
  
  ###########
  phase <- groom.tot$group
  phase <- substr(phase,(nchar(phase)+1)-4,nchar(phase))
  groom.tot$phase <- phase
  
  ## generate a data frame with minutes of grooming (total), elo, # of obs, and group, per individual
  ## phase I
  tmp.minutes <- vector()
  tmp.minutes.given <- vector()
  tmp.minutes.received <- vector()
  tmp.names <- vector()
  tmp.ranks <- vector()
  tmp.obs <- vector()
  tmp.group <- vector()
  for (i in unique(groom.tot$ID1)){
    tmp <- subset(groom.tot, groom.tot$phase == 'SGE1')
    tmp2 <- sum(subset(tmp, tmp$ID1 == i)$min)
    tmp3 <- sum(subset(tmp, tmp$ID1 == i)$ID1.groom.min)
    tmp4 <- sum(subset(tmp, tmp$ID1 == i)$ID2.groom.min)
    tmp.minutes <- cbind(tmp.minutes,tmp2)
    tmp.minutes.given <- cbind(tmp.minutes.given,tmp3)
    tmp.minutes.received <- cbind(tmp.minutes.received,tmp4)
    tmp.names <- cbind(tmp.names,i)
    tmp.ranks <- cbind(tmp.ranks,subset(tmp, tmp$ID1 == i)$ID1.rank[1])
    tmp.obs <- cbind(tmp.obs,subset(tmp, tmp$ID1 == i)$no.obs[1])
    tmp.group <- cbind(tmp.group,subset(tmp, tmp$ID1 == i)$group[1])
  }  
  sge.I.groom.rate.df <- as.data.frame(t(rbind(tmp.names,tmp.ranks,tmp.minutes,tmp.minutes.given,tmp.minutes.received,tmp.obs,tmp.group)))
  names(sge.I.groom.rate.df) <- c('ID','elo','min','min.given','min.received','no.obs','group')
  sge.I.groom.rate.df$elo <- as.numeric(as.character(sge.I.groom.rate.df$elo))
  sge.I.groom.rate.df$min <- as.numeric(as.character(sge.I.groom.rate.df$min))
  sge.I.groom.rate.df$min.given <- as.numeric(as.character(sge.I.groom.rate.df$min.given))
  sge.I.groom.rate.df$min.received <- as.numeric(as.character(sge.I.groom.rate.df$min.received))
  sge.I.groom.rate.df$no.obs <- as.numeric(as.character(sge.I.groom.rate.df$no.obs))
  sge.I.groom.rate.df$phase <- 1
  
  ## now for phase II
  tmp.minutes <- vector()
  tmp.minutes.given <- vector()
  tmp.minutes.received <- vector()
  tmp.names <- vector()
  tmp.ranks <- vector()
  tmp.obs <- vector()
  tmp.group <- vector()
  for (i in unique(groom.tot$ID1)){
    tmp <- subset(groom.tot, groom.tot$phase == 'SGE2')
    tmp2 <- sum(subset(tmp, tmp$ID1 == i)$min)
    tmp3 <- sum(subset(tmp, tmp$ID1 == i)$ID1.groom.min)
    tmp4 <- sum(subset(tmp, tmp$ID1 == i)$ID2.groom.min)
    tmp.minutes <- cbind(tmp.minutes,tmp2)
    tmp.minutes.given <- cbind(tmp.minutes.given,tmp3)
    tmp.minutes.received <- cbind(tmp.minutes.received,tmp4)
    tmp.names <- cbind(tmp.names,i)
    tmp.ranks <- cbind(tmp.ranks,subset(tmp, tmp$ID1 == i)$ID1.rank[1])
    tmp.obs <- cbind(tmp.obs,subset(tmp, tmp$ID1 == i)$no.obs[1])
    tmp.group <- cbind(tmp.group,subset(tmp, tmp$ID1 == i)$group[1])
  }  
  sge.II.groom.rate.df <- as.data.frame(t(rbind(tmp.names,tmp.ranks,tmp.minutes,tmp.minutes.given,tmp.minutes.received,tmp.obs,tmp.group)))
  names(sge.II.groom.rate.df) <- c('ID','elo','min','min.given','min.received','no.obs','group')
  sge.II.groom.rate.df$elo <- as.numeric(as.character(sge.II.groom.rate.df$elo))
  sge.II.groom.rate.df$min <- as.numeric(as.character(sge.II.groom.rate.df$min))
  sge.II.groom.rate.df$min.given <- as.numeric(as.character(sge.II.groom.rate.df$min.given))
  sge.II.groom.rate.df$min.received <- as.numeric(as.character(sge.II.groom.rate.df$min.received))
  sge.II.groom.rate.df$no.obs <- as.numeric(as.character(sge.II.groom.rate.df$no.obs))
  sge.II.groom.rate.df$phase <- 2
  
  sge.groom.rate.df <- rbind(sge.I.groom.rate.df,sge.II.groom.rate.df)
  return(sge.groom.rate.df)
}

for (i in ls(pattern='data_window')){
  tmp <- get(i)
  intros_tmp <- intros
  for (x in names(data)) intros_tmp[which(intros_tmp[,"group"]==x),"num.obs"]=length(levels(as.factor((tmp[[x]][,"obs.no"]))))
  social_tmp <- lapply(tmp,convert.soc)
  groom_tmp <- lapply(social_tmp,function(x){x[which(x[,5]%in%c("GM","G-","ZZ")),]})
  names(groom_tmp) <- names(tmp)
  groom_tmp2 <- groom.function2(groom_tmp)
  agg_tmp <- agg.function(social_tmp)
  assign(value = groom_tmp2, x = paste('groom_df_window_',strsplit(i, split='_')[[1]][3],sep=''))
  assign(value = agg_tmp, x = paste('agg_df_window_',strsplit(i, split='_')[[1]][3],sep=''))
}

test_minute_matrix <- matrix(nrow = 90, ncol = 27)
#Pattern1 <- grep("df_window",names(.GlobalEnv),value=TRUE)
test_list <- list(groom_df_window_1,groom_df_window_2,groom_df_window_3,groom_df_window_4,groom_df_window_5,groom_df_window_6,groom_df_window_7,groom_df_window_8,groom_df_window_9,groom_df_window_10,groom_df_window_11,groom_df_window_12,groom_df_window_13,groom_df_window_14,groom_df_window_15,groom_df_window_16,groom_df_window_17,groom_df_window_18,groom_df_window_19,groom_df_window_20,groom_df_window_21,groom_df_window_22,groom_df_window_23,groom_df_window_24,groom_df_window_25,groom_df_window_26,groom_df_window_27)

for (i in 1:27){
  test_minute_matrix[,i] <- test_list[[i]][,3]
}

test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
  geom_boxplot() +
  xlab('sliding window') +
  ggtitle(label='',subtitle = 'minutes grooming')

for (i in 1:27){
  test_minute_matrix[,i] <- test_list[[i]][,4]
}

test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
  geom_boxplot() +
  xlab('sliding window') +
  ggtitle(label='',subtitle = 'minutes grooming given')

for (i in 1:27){
  test_minute_matrix[,i] <- test_list[[i]][,5]
}

test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
  geom_boxplot() +
  xlab('sliding window') +
  ggtitle(label='',subtitle = 'minutes grooming received')

# agg plots
test_minute_matrix <- matrix(nrow = 90, ncol = 27)

test_list <- list(agg_df_window_1,agg_df_window_2,agg_df_window_3,agg_df_window_4,agg_df_window_5,agg_df_window_6,agg_df_window_7,agg_df_window_8,agg_df_window_9,agg_df_window_10,agg_df_window_11,agg_df_window_12,agg_df_window_13,agg_df_window_14,agg_df_window_15,agg_df_window_16,agg_df_window_17,agg_df_window_18,agg_df_window_19,agg_df_window_20,agg_df_window_21,agg_df_window_22,agg_df_window_23,agg_df_window_24,agg_df_window_25,agg_df_window_26,agg_df_window_27)

for (i in 1:27){
  test_minute_matrix[,i] <- test_list[[i]][,4]
}

test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
  geom_boxplot() +
  xlab('sliding window') +
  ggtitle(label='',subtitle = 'aggs given')

for (i in 1:27){
  test_minute_matrix[,i] <- test_list[[i]][,8]
}

test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
  geom_boxplot() +
  xlab('sliding window') +
  ggtitle(label='',subtitle = 'aggs received')

for (i in 1:27){
  test_minute_matrix[,i] <- test_list[[i]][,6]
}

test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
  geom_boxplot() +
  xlab('sliding window') +
  ggtitle(label='',subtitle = 'contact aggs given')

for (i in 1:27){
  test_minute_matrix[,i] <- test_list[[i]][,7]
}

test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
  geom_boxplot() +
  xlab('sliding window') +
  ggtitle(label='',subtitle = 'noncontact aggs given')

for (i in 1:27){
  test_minute_matrix[,i] <- test_list[[i]][,9]
}

test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
  geom_boxplot() +
  xlab('sliding window') +
  ggtitle(label='',subtitle = 'contact aggs received')

for (i in 1:27){
  test_minute_matrix[,i] <- test_list[[i]][,10]
}

test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
  geom_boxplot() +
  xlab('sliding window') +
  ggtitle(label='',subtitle = 'noncontact aggs received')

Behavioral variable modeling

#### Produce heatmap of FDR sig. genes from all behavioral models on HARDAC 
results_matrix_tmp <- matrix(nrow = 13, ncol = 14)
rownames(results_matrix_tmp) <- c('elo',
                                  'aggDiff',
                                  'Contact agg given',
                                  'Contact agg rec',
                                  'Groom given',
                                  'Groom rec',
                                  'Noncontact agg given',
                                  'Noncontact agg rec',
                                  'Overall agg given',
                                  'Overall agg rec',
                                  'Overall groom',
                                  'PC1.agg',
                                  'PC1.groom')
colnames(results_matrix_tmp) <- c("CD14",
                                  "CD16",
                                  "CD20",
                                  "CD4",
                                  "CD8",
                                  "dexCA",
                                  "dexCA_delta",
                                  "dexGE",
                                  "dexGE_delta",
                                  "gard",
                                  "gard_delta",
                                  "lps",
                                  "lps_delta",
                                  "null")

elo_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_FDR.rds')
elo_CD4 <- elo_flow[[1]]
elo_CD8 <- elo_flow[[2]]
elo_CD14 <- elo_flow[[3]]
elo_CD16 <- elo_flow[[4]]
elo_CD20 <- elo_flow[[5]]
rm(elo_flow)
elo_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_gard_modelRes_fdr')
elo_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/Elo_GARDdelta_modelRes.txt')
elo_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_Elo.txt')
elo_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/Elo_LPSdelta_modelRes.txt')
elo_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_ca_modelRes_fdr')
elo_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_CA_delta_modelRes_fdr')
elo_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_ge_modelRes_fdr')
elo_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_GE_delta_modelRes_fdr')

PC1.agg_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_FDR.rds')
PC1.agg_CD4 <- PC1.agg_flow[[1]]
PC1.agg_CD8 <- PC1.agg_flow[[2]]
PC1.agg_CD14 <- PC1.agg_flow[[3]]
PC1.agg_CD16 <- PC1.agg_flow[[4]]
PC1.agg_CD20 <- PC1.agg_flow[[5]]
rm(PC1.agg_flow)
PC1.agg_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_gard_modelRes_fdr')
PC1.agg_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_GARDdelta_modelRes.txt')
PC1.agg_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_PC1.agg.txt')
PC1.agg_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_LPSdelta_modelRes.txt')
PC1.agg_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_ca_modelRes_fdr')
PC1.agg_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_CA_delta_modelRes_fdr')
PC1.agg_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_ge_modelRes_fdr')
PC1.agg_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_GE_delta_modelRes_fdr')

PC1.groom_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_FDR.rds')
PC1.groom_CD4 <- PC1.groom_flow[[1]]
PC1.groom_CD8 <- PC1.groom_flow[[2]]
PC1.groom_CD14 <- PC1.groom_flow[[3]]
PC1.groom_CD16 <- PC1.groom_flow[[4]]
PC1.groom_CD20 <- PC1.groom_flow[[5]]
rm(PC1.groom_flow)
PC1.groom_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_gard_modelRes_fdr')
PC1.groom_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_GARDdelta_modelRes.txt')
PC1.groom_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_PC1.groom.txt')
PC1.groom_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_LPSdelta_modelRes.txt')
PC1.groom_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_ca_modelRes_fdr')
PC1.groom_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_CA_delta_modelRes_fdr')
PC1.groom_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_ge_modelRes_fdr')
PC1.groom_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_GE_delta_modelRes_fdr')

gc.agg.Diff_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_FDR.rds')
gc.agg.Diff_CD4 <- gc.agg.Diff_flow[[1]]
gc.agg.Diff_CD8 <- gc.agg.Diff_flow[[2]]
gc.agg.Diff_CD14 <- gc.agg.Diff_flow[[3]]
gc.agg.Diff_CD16 <- gc.agg.Diff_flow[[4]]
gc.agg.Diff_CD20 <- gc.agg.Diff_flow[[5]]
rm(gc.agg.Diff_flow)
gc.agg.Diff_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_gard_modelRes_fdr')
gc.agg.Diff_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_GARDdelta_modelRes.txt')
gc.agg.Diff_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.agg.Diff.txt')
gc.agg.Diff_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_LPSdelta_modelRes.txt')
gc.agg.Diff_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_ca_modelRes_fdr')
gc.agg.Diff_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_CA_delta_modelRes_fdr')
gc.agg.Diff_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_ge_modelRes_fdr')
gc.agg.Diff_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_GE_delta_modelRes_fdr')

gc.contact.agg.given_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_FDR.rds')
gc.contact.agg.given_CD4 <- gc.contact.agg.given_flow[[1]]
gc.contact.agg.given_CD8 <- gc.contact.agg.given_flow[[2]]
gc.contact.agg.given_CD14 <- gc.contact.agg.given_flow[[3]]
gc.contact.agg.given_CD16 <- gc.contact.agg.given_flow[[4]]
gc.contact.agg.given_CD20 <- gc.contact.agg.given_flow[[5]]
rm(gc.contact.agg.given_flow)
gc.contact.agg.given_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_gard_modelRes_fdr')
gc.contact.agg.given_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_GARDdelta_modelRes.txt')
gc.contact.agg.given_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.contact.agg.given.txt')
gc.contact.agg.given_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_LPSdelta_modelRes.txt')
gc.contact.agg.given_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_ca_modelRes_fdr')
gc.contact.agg.given_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_CA_delta_modelRes_fdr')
gc.contact.agg.given_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_ge_modelRes_fdr')
gc.contact.agg.given_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_GE_delta_modelRes_fdr')

gc.contact.agg.rec_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_FDR.rds')
gc.contact.agg.rec_CD4 <- gc.contact.agg.rec_flow[[1]]
gc.contact.agg.rec_CD8 <- gc.contact.agg.rec_flow[[2]]
gc.contact.agg.rec_CD14 <- gc.contact.agg.rec_flow[[3]]
gc.contact.agg.rec_CD16 <- gc.contact.agg.rec_flow[[4]]
gc.contact.agg.rec_CD20 <- gc.contact.agg.rec_flow[[5]]
rm(gc.contact.agg.rec_flow)
gc.contact.agg.rec_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_gard_modelRes_fdr')
gc.contact.agg.rec_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_GARDdelta_modelRes.txt')
gc.contact.agg.rec_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.contact.agg.rec.txt')
gc.contact.agg.rec_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_LPSdelta_modelRes.txt')
gc.contact.agg.rec_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_ca_modelRes_fdr')
gc.contact.agg.rec_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_CA_delta_modelRes_fdr')
gc.contact.agg.rec_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_ge_modelRes_fdr')
gc.contact.agg.rec_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_GE_delta_modelRes_fdr')

gc.groom.given_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_FDR.rds')
gc.groom.given_CD4 <- gc.groom.given_flow[[1]]
gc.groom.given_CD8 <- gc.groom.given_flow[[2]]
gc.groom.given_CD14 <- gc.groom.given_flow[[3]]
gc.groom.given_CD16 <- gc.groom.given_flow[[4]]
gc.groom.given_CD20 <- gc.groom.given_flow[[5]]
rm(gc.groom.given_flow)
gc.groom.given_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_gard_modelRes_fdr')
gc.groom.given_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_GARDdelta_modelRes.txt')
gc.groom.given_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.groom.given.txt')
gc.groom.given_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_LPSdelta_modelRes.txt')
gc.groom.given_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_ca_modelRes_fdr')
gc.groom.given_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_CA_delta_modelRes_fdr')
gc.groom.given_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_ge_modelRes_fdr')
gc.groom.given_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_GE_delta_modelRes_fdr')

gc.groom.rec_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_FDR.rds')
gc.groom.rec_CD4 <- gc.groom.rec_flow[[1]]
gc.groom.rec_CD8 <- gc.groom.rec_flow[[2]]
gc.groom.rec_CD14 <- gc.groom.rec_flow[[3]]
gc.groom.rec_CD16 <- gc.groom.rec_flow[[4]]
gc.groom.rec_CD20 <- gc.groom.rec_flow[[5]]
rm(gc.groom.rec_flow)
gc.groom.rec_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_gard_modelRes_fdr')
gc.groom.rec_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_GARDdelta_modelRes.txt')
gc.groom.rec_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.groom.rec.txt')
gc.groom.rec_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_LPSdelta_modelRes.txt')
gc.groom.rec_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_ca_modelRes_fdr')
gc.groom.rec_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_CA_delta_modelRes_fdr')
gc.groom.rec_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_ge_modelRes_fdr')
gc.groom.rec_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_GE_delta_modelRes_fdr')

gc.noncontact.agg.given_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_FDR.rds')
gc.noncontact.agg.given_CD4 <- gc.noncontact.agg.given_flow[[1]]
gc.noncontact.agg.given_CD8 <- gc.noncontact.agg.given_flow[[2]]
gc.noncontact.agg.given_CD14 <- gc.noncontact.agg.given_flow[[3]]
gc.noncontact.agg.given_CD16 <- gc.noncontact.agg.given_flow[[4]]
gc.noncontact.agg.given_CD20 <- gc.noncontact.agg.given_flow[[5]]
rm(gc.noncontact.agg.given_flow)
gc.noncontact.agg.given_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_gard_modelRes_fdr')
gc.noncontact.agg.given_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_GARDdelta_modelRes.txt')
gc.noncontact.agg.given_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.noncontact.agg.given.txt')
gc.noncontact.agg.given_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_LPSdelta_modelRes.txt')
gc.noncontact.agg.given_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_ca_modelRes_fdr')
gc.noncontact.agg.given_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_CA_delta_modelRes_fdr')
gc.noncontact.agg.given_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_ge_modelRes_fdr')
gc.noncontact.agg.given_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_GE_delta_modelRes_fdr')

gc.noncontact.agg.rec_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_FDR.rds')
gc.noncontact.agg.rec_CD4 <- gc.noncontact.agg.rec_flow[[1]]
gc.noncontact.agg.rec_CD8 <- gc.noncontact.agg.rec_flow[[2]]
gc.noncontact.agg.rec_CD14 <- gc.noncontact.agg.rec_flow[[3]]
gc.noncontact.agg.rec_CD16 <- gc.noncontact.agg.rec_flow[[4]]
gc.noncontact.agg.rec_CD20 <- gc.noncontact.agg.rec_flow[[5]]
rm(gc.noncontact.agg.rec_flow)
gc.noncontact.agg.rec_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_gard_modelRes_fdr')
gc.noncontact.agg.rec_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_GARDdelta_modelRes.txt')
gc.noncontact.agg.rec_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.noncontact.agg.rec.txt')
gc.noncontact.agg.rec_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_LPSdelta_modelRes.txt')
gc.noncontact.agg.rec_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_ca_modelRes_fdr')
gc.noncontact.agg.rec_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_CA_delta_modelRes_fdr')
gc.noncontact.agg.rec_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_ge_modelRes_fdr')
gc.noncontact.agg.rec_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_GE_delta_modelRes_fdr')

gc.overall.agg.rec_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_FDR.rds')
gc.overall.agg.rec_CD4 <- gc.overall.agg.rec_flow[[1]]
gc.overall.agg.rec_CD8 <- gc.overall.agg.rec_flow[[2]]
gc.overall.agg.rec_CD14 <- gc.overall.agg.rec_flow[[3]]
gc.overall.agg.rec_CD16 <- gc.overall.agg.rec_flow[[4]]
gc.overall.agg.rec_CD20 <- gc.overall.agg.rec_flow[[5]]
rm(gc.overall.agg.rec_flow)
gc.overall.agg.rec_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_gard_modelRes_fdr')
gc.overall.agg.rec_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_GARDdelta_modelRes.txt')
gc.overall.agg.rec_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.overall.agg.rec.txt')
gc.overall.agg.rec_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_LPSdelta_modelRes.txt')
gc.overall.agg.rec_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_ca_modelRes_fdr')
gc.overall.agg.rec_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_CA_delta_modelRes_fdr')
gc.overall.agg.rec_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_ge_modelRes_fdr')
gc.overall.agg.rec_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_GE_delta_modelRes_fdr')

gc.overall.agg.given_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_FDR.rds')
gc.overall.agg.given_CD4 <- gc.overall.agg.given_flow[[1]]
gc.overall.agg.given_CD8 <- gc.overall.agg.given_flow[[2]]
gc.overall.agg.given_CD14 <- gc.overall.agg.given_flow[[3]]
gc.overall.agg.given_CD16 <- gc.overall.agg.given_flow[[4]]
gc.overall.agg.given_CD20 <- gc.overall.agg.given_flow[[5]]
rm(gc.overall.agg.given_flow)
gc.overall.agg.given_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_gard_modelRes_fdr')
gc.overall.agg.given_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_GARDdelta_modelRes.txt')
gc.overall.agg.given_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.overall.agg.given.txt')
gc.overall.agg.given_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_LPSdelta_modelRes.txt')
gc.overall.agg.given_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_ca_modelRes_fdr')
gc.overall.agg.given_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_CA_delta_modelRes_fdr')
gc.overall.agg.given_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_ge_modelRes_fdr')
gc.overall.agg.given_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_GE_delta_modelRes_fdr')

gc.overall.groom_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_FDR.rds')
gc.overall.groom_CD4 <- gc.overall.groom_flow[[1]]
gc.overall.groom_CD8 <- gc.overall.groom_flow[[2]]
gc.overall.groom_CD14 <- gc.overall.groom_flow[[3]]
gc.overall.groom_CD16 <- gc.overall.groom_flow[[4]]
gc.overall.groom_CD20 <- gc.overall.groom_flow[[5]]
rm(gc.overall.groom_flow)
gc.overall.groom_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_gard_modelRes_fdr')
gc.overall.groom_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_GARDdelta_modelRes.txt')
gc.overall.groom_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.overall.groom.txt')
gc.overall.groom_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_LPSdelta_modelRes.txt')
gc.overall.groom_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_ca_modelRes_fdr')
gc.overall.groom_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_CA_delta_modelRes_fdr')
gc.overall.groom_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_ge_modelRes_fdr')
gc.overall.groom_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_GE_delta_modelRes_fdr')


# template code for inputting # sig genes
#countFlag <- 1
#for (i in ls(pattern = 'elo_')){
#  tmp <- get(i)
#  results_matrix_tmp[1,countFlag] <- sum(tmp[,ncol(tmp)] < 0.1)
#  countFlag <- countFlag + 1
#}
#results_matrix_tmp[1,14] <- sum(elo_lps$fdr_Elo_at_NC < 0.1)


behaviorNames <- c("elo_",
                   "gc.agg.Diff_",
                   "gc.contact.agg.given_",
                   "gc.contact.agg.rec_",
                   "gc.groom.given_",
                   "gc.groom.rec_",
                   "gc.noncontact.agg.given_",
                   "gc.noncontact.agg.rec_",
                   "gc.overall.agg.given_",
                   "gc.overall.agg.rec_",
                   "gc.overall.groom_",
                   "PC1.agg_",
                   "PC1.groom_")

countFlag_beh <- 1
countFlag_dataset <- 1
for (b in behaviorNames){
  tmp_filelist <- ls(pattern=b)
  countFlag_dataset <- 1
  for (i in tmp_filelist){
    tmp <- get(i)
    results_matrix_tmp[countFlag_beh,countFlag_dataset] <- sum(tmp[,ncol(tmp)] < 0.1)
    countFlag_dataset <- countFlag_dataset + 1
  }
  countFlag_beh <- countFlag_beh + 1
}

countFlag_beh <- 1
for (i in behaviorNames){
  tmp <- paste0(i, 'lps', sep='')
  tmp <- get(tmp)
  results_matrix_tmp[countFlag_beh,14] <- sum(tmp[,7] < 0.1)
  countFlag_beh <- countFlag_beh + 1
}

results_matrix <- results_matrix_tmp[2:13,]
for (i in 1:ncol(results_matrix)){
  results_matrix[,i] <- log2((results_matrix[,i] +1) / (results_matrix_tmp[1,i] +1))
}

colnames(results_matrix) <- c("CD14 flow-sorted (0)",
                              "CD16 flow-sorted (1635)",
                              "CD20 flow-sorted (83)",
                              "CD4 flow-sorted (326)",
                              "CD8 flow-sorted (23)",
                              "Dex CA (1)",
                              "Dex CA_delta (1)",
                              "Dex GE (1)",
                              "Dex GE_delta (1)",
                              "Gard (3517)",
                              "Gard_delta (0)",
                              "LPS (6234)",
                              "LPS_delta (4957)",
                              "NULL (4487)")


pheatmap(as.matrix(t(results_matrix)), color=colorRampPalette(c("lightblue4", "white", 'mediumpurple4'))(100), border_color = 'white', scale = 'row', legend = T, treeheight_row = 0, treeheight_col = 0, main = 'log2(FC) # of sig. genes over Elo', cluster_rows = 0, cluster_cols = 1)